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Roundoff errors cannot be avoided when implementing numerical programs with finite precision. The ability 
to reason about rounding is especially important if one wants to explore a range of potential representations, 
for instance for FPGAs or custom hardware implementations. This problem becomes challenging when the 
program does not employ solely linear operations, and non-linearities are inherent to many interesting 
computational problems in real-world applications. 

Existing solutions to reasoning possibly lead to either inaccurate bounds or high analysis time in the 
presence of nonlinear correlations between variables. Furthermore, while it is easy to implement a straight¬ 
forward method such as interval arithmetic, sophisticated techniques are less straightforward to implement 
in a formal setting. Thus there is a need for methods which output certificates that can be formally validated 
inside a proof assistant. 

We present a framework to provide upper bounds on absolute roundoff errors of fioating-point nonlinear 
programs. This framework is based on optimization techniques employing semidefinite programming and 
sums of squares certificates, which can be checked inside the Coq theorem prover to provide formal roundoff 
error bounds for polynomial programs. Our tool covers a wide range of nonlinear programs, including 
polynomials and transcendental operations as well as conditional statements. We illustrate the efficiency 
and precision of this tool on non-trivial programs coming from biology, optimization and space control. Our 
tool produces more accurate error bounds for 23 % of all programs and yields better performance in 66 % 
of all programs. 
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1. INTRODUCTION 

Constructing numerical programs which perform accurate computation turns out to be dif¬ 
ficult, due to finite numerical precision of implementations such as fioating-point or fixed- 
point representations. Finite-precision numbers induce roundoff errors, and knowledge of 
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the range of these roundoff errors is required to fulfill safety criteria of critical programs, as 
typically arising in modern embedded systems such as aircraft controllers. Such a knowledge 
can be used in general for developing accurate numerical software, but is also particularly 
relevant when considering migration of algorithms onto hardware (e.g. FPGAs). The ad¬ 
vantage of architectures based on FPGAs is that they allow more flexible choices in number 
representations, rather than limiting the choice between IEEE standard single or double 
precision. Indeed, in this case, we benefit from a more flexible number representation while 
still ensuring guaranteed bounds on the program output. 

To obtain lower bou nds on roundoff err ors, one can rely on testing approaches, su ch as 

or under-approximation tools (e.g. s3fp 


meta-heuristic search 
Here 


Borges et al. 2012 


Ghi- 


ang et al. 2014 ). Here, we are interested in efficiently handling the complementary over¬ 
approximation problem, namely to obtain precise upper bounds on the error. This problem 
boils down to finding tight abstractions of linearities or non-linearities while being able 
to bound the resulting approximations in an efficient way. For computer programs con¬ 
sisting of linear operations, automatic error analysis can be obtained with well-studied 
optimi zation techniques ba sed on SAT/SMT solvers [Haller et al. 2012| and affine arith¬ 


metic Delmas et al. 2009 . However, non-linear operations are key to many interesting 


computational problems arising in physics, biology, controller implementations and global 
optimization. Recently, two promising frameworks have been designed to provide upper 
bounds for roundoff errors of nonlinear programs. The corresponding algorithms rely on 


implemented in the FP Taylor tool, and on 
implemented in the 


Taylor-interval methods Solovyev et al. 2015 

combining SMT with interval arithmetic Darulova and Kuncak 2014 
Rosa real compiler. 

The complexity of the mathematics underlying techniques for nonlinear reasoning, and 
the intricacies associated with constructing an efficient implementation, are such that a 
means for independent formal validation of results is particularly desirable. The Rosa tool is 
based on theoretical results that should provide sound over-approximations of error bounds. 
While Rosa relies on an SMT solver capable of generating unsatisfiability proof witnesses 
(thus allowing independent soundness checking), it does not formally verify these certifi¬ 
cates inside a proof assistant. To the best of our knowledge, FPTaylor and Gappa are the 
only academic software tools that can pro duce formal proof certifica tes. For FPTaylor, this 
is based on the framewor k developed in Solovyev and Hales 2013 to verify nonlinear in¬ 


equalities in Hol-light Harrison 1996 using Taylor-interval methods. However, most of 


computation performed in the informal optimization procedure ends up being redone inside 
the Hol-light proof assistant, yielding a formal verification which may be computationally 
demanding. 

The aim of this work is to provide a formal framework to perform automated precision 
analysis of computer programs that manipulate finite-precision data using nonlinear opera¬ 
tors. For such programs, guarantees can be provided with certified programming techniques. 
Semidefinite programming (SDP) is relevant to a wide range of mathematical fields, includ¬ 
ing combinatorial optimization, control t heory and mat rix completion. In 2001, Lasserre 
introduced a hierarchy of SDP relaxations Lasserre 2001 for approximating polynomial in- 
fima. Our method to b ound the error i s a decision procedure based on a specialized variant of 
the Lasserre hierarchy Lasserre 2006| . The procedure relies on SDP to provide sparse sum- 
of-squares decompositions of nonnegative polynomials. Our framework handles polynomial 
program analysis (involving the operations -F, x, —) as well as extensions to the more general 
class of semialgebraic and transcendental programs (involving min, max, arctan, exp), 
following the approximation scheme described in [Magron et al. 2015a . 
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1.1. Overview of our Method 

We present an overview of our method and of the capabilities of related techniques, using 
an example. Consider a program implementing the following polynomial expression /: 


/(x) := a;2 X X5 + X3 X ^6 - X2 X X3 - ^5 X xq 
+XI X (—X1 + X2 + X:i — X4^ + X^ + Xq) , 

where the six-variable vector x := {xi,X2,X3,X4,X5,Xq) is the input of the program. For 
this example, assume that the set X of possible input values is a product of closed intervals: 
X = [4.00,6.36]®. This function / together with the set X a ppear in many inequalities 
arising from the the proof of the Kepler Conjecture Hales 2006 , yielding challenging global 
optimization problems. 

The polynomial expression / is obtained by performing 15 basic operations (1 negation, 
3 subtractions, 6 additions and 5 multiplications). When executing this program with a set 
of floating-point numbers x := {xi,X 2 ,x^, X 4 _, £ 5 , xe) G X, one actually computes a floating¬ 
point result /, where all operations -|-,—,x are replaced by the respectively associated 
floating-point opera tions 0,0,0 . The results of these operations comp ly with IEEE 754 
standard arithmetic IEEE 2008] (see relevant background in Section |2.1|). Here, for the sake 
of clarity, we do not consider real input variables but we do it later on while performing 
detailed comparison (see Section]^. For instance, (in the absence of underflow) one can 
write X 2 ®x^ = {x 2 x a; 5 )(l-|-ei), by introducing an error variable ei such that —e < ei < e, 
where the bound e is the machine precision (e.g. e = 2“^^ for single precision). One would 
like to bound the absolute roundoff error |r(x,e)| := |/(x,e) — /(x)| over all possible input 
variables x G X and error variable ei,...,ei 5 G [—e, ej. Let us deflne E := [—e, and 
K := X X E. Then our bound problem can be cast as finding the maximum r* of | r | over 
K, yielding the following nonlinear optimization problem: 


r* := max |r(x,e)| 

(x,e)6K 


= max{— min r(x,e), max r(x,e)} , 

(x,e)eK (x,e)eK 


( 1 ) 


One can directly try to solve th ese two polynomial optimiza tion problems using classical 

one can also decompose the 
and a nonlinear term 


SDP relaxations Lasserre 2001 
error term r as t 


As in 

e sum of a term l(x,e 


[Solovye 


V et al. 2015 
which is affine w.r.t. e 
/i(x,e) := r(x,e) — Z(x,e). Then the triangular inequality yields: 


< max U(x,e)| 

(x,e)GK 


max |/i(x,e)| 

(x.e)GK' 


( 2 ) 


It follows for this example that /(x,e) = X2X5ei + x^XQe2 + {x2X^ + x^XQ)e^ + - ■ •-|-/(x)ei5 = 
^15 „ ^ 2:33:6, Si5(x) := /(x). The Symbolic Taylor 

consists of using a simple branch and bound 


Si(x)ei, with Si 
Expansions method 


X) := X2X5,S2(X} 


Solovyev et al. 2015 

algorithm based on interval arithmetic to compute a rigorous interval enclosure of each 
polynomial Si, i = 1,..., 15, over X and Anally obtain an upper bound of jZj -I- \h\ over K. 
In contrast, o ur method use s sparse semideflnite relaxations for polynomi al optimization 
(deriv ed from [Lasserre 2006 ) to bound I and basic interval arithmetic as in Solovyev et al. 
2015 to bound \h\ (i.e. we use interval arithmetic to bound second-order error terms in the 
multivariate Taylor expansion of r w.r.t. e). 

The following comparison results have been obtained on an Intel Core i7-5600U CPU 
(2.60 GHz). All execution times have been computed by averaging over five runs. 

— A direct attempt to solve the two po lynomial problems occ urring in Equation Q fails 
as the SDP solver (in our case Sdpa Yamashita et al. 2010 ) runs out of memory. 
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Using our method implemented in the Real2Float tool, one obtains an upper bound of 
760e for |/| + |/i| over K in 0.15 seconds. This bound is provided together with a certificate 
which can be formally checked inside the COQ proof assistant in 0.20 seconds. 

After normalizing the polynomial expression and using basic interval arithmetic, one 
obtains 8 times more quickly a coarser bound of 922e. 

Symbolic Taylor expansions implemented in FPTaylor 
more precise bound of 721e, but the analysis time is 


Solovyev et al. 2015 provide a 
times slower than with our 


28 

implementation. Formal verification of this bound inside the Hol-light proof assistant 
takes 27.7 seconds, which is 139 times slower than proof checking with Real2Float 
inside COQ. One can obtain an even more precise bound of 528e (but 37 times slower 
than with our implementation) by turning on the improved rounding model of FPTaylor 
and limiting the number of branch and bound iterations to 10000. The drawback of this 
bound is that it cannot be formally verified. 


— Finally, a sligthly coarser bound of 762e is obtained with the Rosa real compiler Darulova 


and Kuncak 2014 , but the analysis is 19 times slower than with our implementation and 
we cannot get formal verification of this bound. 

1.2. Related Works 

SMT solvers allow analysis of programs with various semantics or specifications but are 
limited for th e manipulation of problems involving nonlinear arithmetic. Several solvers, 
including Z3 De Mour a and Bjprner 2008 , provide partial support for the IEEE floating¬ 
point standard” Rummer and Wahl 2010 


They suffer fr om a lack of scalability when used 
for roundoff error analysis in isolation (as emphasized in [Darulova and Kuncak 2014|) , but 
can be integrated into existing frameworks, e.g. FPhile Faganeili and Ahrendt 2013 . The 


procedure in Gao et al. 2013 


can solve SMT problems over the real numbers, using interval 


constraint propa gation, buFhas not yet been applied to quantification of roundoff error. 


The Rosa tool Darulova and Kuncak 2014 provides a way to compile functional SCALA 


programs involving semialgebraic functions and conditional statements. The tool uses affine 
arithmetic to provide sound over-approximations of roundoff errors, allowing for generation 
of finite precision implementations which fulfill the required precision given as input by 
the user. This tool thus relies on abstract interpration but bounds of the affine expressions 
are provided through an optimization procedure based on SMT. In our case, we use the 
same rounding model but provide approximations which are affine w.r.t. the additional 
error variables and nonlinear w.r.t. the input variables. Instead of using SMT, we bound 
the resulting expressions with optimization techniq ues based on semidefinite programming. 
Abstract interpretation Cousot and Cousot 1977] has been extensively used in the con¬ 


text of static analysis to provide sound over-approximations, called abstractions, of the 
sets of values taken by program variables. The effects of variable assignments, guards and 
conditional branching statements are handled with several domain specifi c operators (e .g. 
inclusion, meet and join). Well studied abstract domains include intervals Moore 1962 as 
well a s more co mplicated fr ameworks b ased on affine arithm etic [Stolfl and de F'igueiredo 
2003 , octogons Mine 200*^, zonotop es Ghorbal et al. 2010 , polyhedra Ghen et al. 2008 , 
interval polyhedra [Ghen et al. 2009| , some of them being irnplemented inside a tool called 
Apron [Jeannet and Mine 


im 


Abstract domains provide sound over-approximations 
of program e xpressions, and allow upper bo unds on roundoff error to be computed. The 
Gappa tool Daumas and Melquiond 2010 relies on interval abstract domains with an 


extension to athne domains [Linderman et al. 2010j, to reason about roundoff errors. As 


demonstrated in Solovyev et al. 2015 , the bounds obtained inside Gappa are often coarser 
than other methods. Formal guarantees can be provided as Gappa benefits from an inter¬ 


face with COQ while making use of interval libraries Melquiond 2012 relying on formalized 
floating-points Boldo and Melquiond 2011 . The static analysis c ommercial tool Fluct uat 
(with a free academic version) relies on affine abstract domains Blanchet et al. 2003 and 
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Table I. Comparison of roundoff error tools w.r.t. expressiveness. 


Feature 

Real2Float 

Rosa 

FPTaylor 

Gappa 

Fluctuat 

Basic FP operations/formats 

7 


y 

y 

y 

Special values (iboo, MaM) 




y 

y 

Improved rounding model 



y 

y 

y 

Input uncertainties 

y 

y 

y 

y 

y 

Transcendental functions 

y 


y 



Discontinuity errors 

y 

y 



y 

Proof certificates 

y 


y 

y 



techniques which are very similar to the ones in Rosa, including interval subdivision. This 
tool does not perform optimization but uses forward computation to analyze floating-point 
prog rams written in C. F urthermore, Fluctuat also has a procedure for discontinuity er¬ 
rors [Ghorbal et al. 20l0 . The Gappa and Fluctuat tools use a different rounding model 
(also available as an option inside FPTaylor) based on a piecewise constant absolute error 
bound. This is more precise than the simple rounding model used in our framework but 
requires (possibly) extensive use of a branch and bound algorithm as each interval has to 
be su bdivided in intervals [2", 2"+^] for several values of the integer n. In Solovyev et al. 

with 


2015 , the authors provide a table (Table 1) comparing relevant features of FPTaylor 
three other tools (Rosa, Gappa and Fluctuat), performing roundoff error estimation. In 
a similar fashion, we summarize the main features related to our tool Real2Float and the 
same four above-mentioned tools used for our further benchmark comparisons w.r.t. their 
expressiveness in Table |T] 

Gomputing sound bounds of nonlinear expressions is mandatory to perform formal analy¬ 
sis of finite precision implementations and can be performed with various optimization tools. 
In the polynomial case, alternative approaches to semidefinite relaxations are based on de¬ 
composition in the multivariate Bernstein basis. Formal verification o f bounds obtained with 
this d ecomposition has been investigated by Munoz and Narkawicz Munoz and Narkawicz 
2013| in the PVS theorem prover. We are not aware of any work based on these techniques 
whicn can quantify roundoff errors. Another decomposition of nonnegative polynomials into 


SOS certificates consists in using the Krivine-Handelman Krivine 1964 Handelman 1988 


representation and boils down to solving linear programming (LP) relaxations. In our case) 
we use a different representation, leading to solve SDP relaxations. The Krivine-Handelman 


representation has been used in Boland and Gonstantinides 2010 to compute roundoff er¬ 


ror bounds. LP relaxations ofteri provide coarser bounds than SDP relaxations and it has 


been proven in Lasserre 2009 that generically finite converg ence does no t occur for convex 
problems, with the exception of the linear case. The work in |Roux 2015] focuses on formal¬ 


ization of roundoff errors bounds related to positive definiteness verification. Branch and 


bound methods with Taylor models Berz and Makino 2009 are not restricted to polynomial 
systems and have been formalized Solovyev and Hales 2013| to solve nonlinear inequalities 
occur ring in the proof of Kepler Conjecture. Symbolic 'i'aylor Expansions Solovyev et al. 
2015 have been implemented in the FPTaylor tool to compute formal bounds of roundoff 


errors for programs involving both polynomial and transcendental functions. 


1.3. Contributions 

Our key contributions can be summarized as follows: 

— We present an optimization algorithm providing sound over-approximations for roundoff 
errors of floating-poin t nonlinear pro grams. This algorithm is based on sparse sums of 

In comparison with other methods, our algorithm 


Lasserre 2006 


squares programming 

allows us to obtain tighter upper bounds, while overcoming scalability and numerical is¬ 


sues inherent in SDP solvers [Todd 20M . Our algorithm can currently handle programs 
implementing polynomial functions, but also involving non-polynomial components, in- 
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eluding either semialgebraic or transcendental operations (e.g. /, ^^arctan, exp), as well 
as conditional statements. Programs containing iterative or while loops are not currently 
supported. 

Our framework is fully implemented in the Real2Float tool. Among several features, the 
tool can optionally perform formal verifi cation of ro undoff error bounds for polynomial 
programs, inside the COQ proo f assistant Co g 2016 . The most recent software release of 
Real2Float provides OCaml OCaml 2015 and COQ libraries and is freely av ailable F 
Our i mplementation tool is built on top of the NLCertify verification system Magron 
20141. Precision and efficiency of the tool are evaluated on several benchmarks coming 


from the existing literature. Numerical experiments demonstrate that our method com¬ 


petes well w ith recent approaches relying on Taylor-interval ap proximations Solov^w 
et al. 2015 or combining SMT solvers with affine arithmetic Darulova and Kuncak 


2014 . Wealso compared o ur tool with Gappa [Daumas and Melquiond 2010 and Fluc- 


TUAT 


Delmas et al. 2009 


The paper is organized as follows. In Section we present mandatory background on 
roundoff errors du e to finite precision arithmetic before describing our nonlinear program 
semantics (Section 2.11. Then we recall how to p erform certified polynomial optimization 
based on semidefinite programming (Section |2.2[ ) and how to obt ain formal bounds while 
checking the certificates inside the COQ proof assistant (Section 2.31. Section contains 
the main contribution of the paper, namely how to compute tight over-approximations for 
roundoff errors of nonlinear programs with sparse semidefinite relaxations. Finally, Section]^ 
is devoted to the evaluation of our nonlinear verification tool Real2Float on benchmarks 
arising from control systems, optimization, physics and biology, as well as comparisons with 
the tools FPTaylor, Rosa, Gappa and Fluctuat. 


2. PRELIMINARIES 

2.1. Program Semantics and Floating-point Numbers 

We support conditional code without procedure calls or loops. Despite these restrictions, we 
can consider a wide range of nonlinear programs while assuming that important numerical 
calculations can be expressed in a loop-free manner. Our programs are encoded in an ML- 
like language: 

let box_prog xi ... a:„ = [(oi, 6i);...; (a„, 6„)]; ; 
let obj_prog .. . x„ = [(/(x), eReai 2 Fioat)] ; ; 

let cstr_prog xi .. . = [ 5 i(x);...; ^^(x)] ; ; 

let uncert_prog xi .. . = [rti;...; u„] ; ; 

Here, the first line encodes interval floating-point bound constraints for input variables, 

namely x := (xi,...,x„) G [ai,6i] x ••• x The second line provides the function 

/(x) as well as the total roundoff error bound eRaai 2 Fioat- Then, one encodes polynomial 
nonnegativity constraints over the input variables, namely gi(x) > 0, ...,( 7 fc(x) > 0. Fi¬ 
nally, the last line allows the user to specify a numerical constant u, to associate a given 
uncertainty to the variable Xi, for each i = 1,... ,n. 

The type of numerical constants is denoted by C. In our current implementation, the 
user can choose either 64 bit floating-point or arbitrary-size rational numbers. This type 
C is used for the terms eReai 2 Fioat, mi, ■ ■ •, Un, oi,..., a„, 5i,..., The inductive type of 
polynomial expressions with coefficients in C is pExprC defined as follows: 

type pexprC = Pc of C I Px of positive 
I Psub of pexprC * pexprC I Pneg of pexprC 
I Padd of pexprC * pexprC 


^ forge.ocamlcore.org/frs/?group_id=351 
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I Pmul of pexprC * pexprC 

The constructor Px takes a positive integer as argument to represent either an input or local 
variable. The inductive type nlexpr of nonlinear expressions (such as /(x)) is defined as 
follows: 


type nlexpr = 

I Pol of pexprC I Neg of nlexpr 
I Add of nlexpr * nlexpr 
I Mul of nlexpr * nlexpr 
I Sub of nlexpr * nlexpr 

I Div of nlexpr * nlexpr I Sqrt of nlexpr 
I Transc of transc*nlexpr 
I IfThenElse of pexprC*nlexpr*nlexpr 
I Let of positive * nlexpr * nlexpr 


The type transc corresponds to a dictionary D of special functions. In our case T> := 
{exp,log,cos,sin,tan,arccos,arcsin,arctan}. For instance, the term Transc (exp, /(x)) 
represents the program implementing exp(/(x)). 

Given a polynomial expression p and two nonlinear expressions / and g, the term 
IfThenElse (p(x), /(x), g(x)) represents the conditional program implementing the ex- 
pressiorj^ if (p(x) > 0) /(x) else ^(x). The constructor Let allows us to define local 
variables in an ML fashion, e.g. let ti = 331.4 + 0.6 * T in —ti * v/{{ti + u) * (ti + u)) 
(part of the dopplerl program considered in Section]^. 

Finally, one obtains rounded nonlinear expressions using a recursive procedure round, 
defined according to Equation ^ and Equation ffl. Rounded expressions are supported 
inside conditions. When an uncertainty Ui is specined for an input variable Xi, the corre¬ 
sponding rounded expression is given by Xi (1 -T e), with | e | < iti, the uncertainty Ui being 
a relative error. _ 

We adopt the standard practice Higham 2002] to approximate a real number x with its 
closest floating-point representation x = a::(l-Te), with |e| is less than the machine precision 
e. In the sequel, we neglect both overflow and denormal range values. The operator ' is called 
the rounding operator and can be selected among rounding to nearest, rounding toward zero 
(resp. ±oo). In the sequel, we assume rounding to nearest. The scientific notation of a binary 
(resp. decimal) floating-point number i is a triple (s, sig, exp) consisting of a sign bit s, a 
significand sig G [1,2) (resp. [1,10)) and an exponent exp, yielding numerical evaluation 
{-!)“ sig2^^P (resp. (-1)® sz^ 10 ®"^p). 

The value of e actually gives the upper bound on the relative floating-point error and is 
equal to where prec is called the precision, referring to the number of significand bits 

used. For single precision floating-point, one has prec = 24. For double (resp. quadruple) 
precision, one has prec = 53 (resp. prec = 113). Let K denote the set of real numbers and 
F the set of binary floating-point numbers. 

For each real-valued operation bopjj G {-T,—,x,/}, the result of the corresponding 
floating-point operation bopj. G |0, Q, 0, 0| satisfies the following when complying with 

IEEE 2008] (without overflow, underflow and denormal 


IEEE 754 standard arithmetic 
occurrences): 


bopj. (x, y) = bopjj {x, y) (1 -T e) 


|< e = 2“ 


( 3 ) 


Other operations include special functions taken from P, containing the unary functions tan, 
arctan, cos, arccos, sin, arcsin, exp, log, (•)^ with r € E\{0}. For the corresponding 


^Our general framework could theoretically handle nested if statements. However, our current implementa¬ 
tion is limited to programs involving single top level conditional statements 
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floating-point evaluation satisfles 


frix) = fR{x){l + e) , 


I e |< eifw) . 


( 4 ) 


The value of the relative error bound e(/R) differs from the machine precision e in Equa¬ 
tion © a nd has to be properly adjusted on a per-operator basis. We refer the interested 
reader to [Bingham and Leslie-Hurd 2014| for relative error bound verification of transcen¬ 


dental functions (see also Harrison 2000 for formalization in Hol-light). 


2.2. SDP Relaxations for Polynomial Optimization 

The sums of squares method involves approximation of polynomial inequality constraints 
by sums of squares (SOS) equality constraints. Here we recall mandatory background about 
SOS. We apply this method in Section]^ to solve the problems of Equation Q when the 
nonlinear function r is a polynomial. In the sequel, let us denote by n the number of initial 
variables of the polynomial optimization problem and by k the number of optimization 
constraints. 


2.2.1. Sums of squares certificates and SDP. First we recall basic facts about generation of 
SOS certificates for polynomial optimization, using semideflnite programming, which can 


be found in texts such as 
and by K 2 d[x] the restriction of. 
set of SOS polynomials: 


Lasserre 2001 . Denote by K[x] the vector space of polynomials 


to polynomials of degree at most 2d. Let us define the 


withg, gK[x]| , 


( 5 ) 


as well as its restriction S 2 d[x] := S[x] p|IR 2 d[x] to polynomials of degree at most 2d. For 
instance, the following bivariate polynomial (t(x) := 1 -I- (xf — lies in S 4 [x] C K 4 [x]. 

Optimization methods based on SOS use the implication r G E[x] Vx G K", r(x) > 
0, i.e. the inclusion of E[x] in the set of nonnegative polynomials. 

The underlying reason for usin g SOS polynomials is that optimizing over positive polyno¬ 
mials is NP Hard Laurent 2009 . Thus, one would like to replace such positivity constraints 


by more tractable ones, and in particular the SOS decompositions admitted by positive 
polynomials provide a suitable alternative: when fixing the degree of such decompositions, 
the resulting relaxed problem becomes more tractable. 

Given r G K[x], one considers the following polynomial minimization problem: 


:= inf I r(x) 

xeB" 


X G K} 


( 6 ) 


where the set of constraints K C K" is defined by 

K := {x G M" : 5i(x) > 0,... ,5ffe(x) > 0} , 

for polynomial functions gi,... ,gk- The set K is called a basic semialgebraic set. Member¬ 
ship of semialgebraic sets is ensured by satisfying conjunctions of polynomial nonnegativity 
constraints. 


Remark 2.1. When the input variables satisfy interval constraints x G [ai,5i] x 
• • • X [a„,6„] then one can easily show that there exists some integer M > 0 such that 
M — X]r=i sequel, we assume that this nonnegativity constraint appears 

explicitly in the definition of K. Such an assumption is mandatory to prove the convergence 
of semideflnite relaxations recalled in Theorem 12.31 

In general, the objective function r and the set of constraints K can be nonconvex, which 
makes Problem ([^ difficult to solve in practice. One can rewrite Problem ^ as the equiv- 
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alent maximization problem: 

r* •= sup{ /r : r(x) — /i > 0 , Vx S K } . (7) 

/^GR 

Now we outline how to handle the nonnegativity constraint r — ^ > 0. Given a nonnegative 
polynomial p G K[x], the existence of an SOS decomposition p = valid over K”, is 

equivalent to the existence of a symmetric real matrix Q, a solution of the following linear 
matrix feasibility problem: 

r(x) = md(x)T Qmd(x), Vx e K”, (8) 


where md(x) := ... ,Xn,Xi,xiX 2 , ■ ■ ■ and the matrix Q has only nonnegative 

eigenvalues. Such a matrix Q is called positive semidefinite. The vector (resp. matrix 
Q) has a size (resp. dimension) equal to := Problem ([^ can be handled with 


semidefinite programming (SDP) solvers, such as Mosek Andersen and Andersen 2000 


or 


Yamashita et al. 2010| (see [Vandenberghe and Poyd 1994 tor specific background 


SDPA 

about SDP). Then, one computes the “LDL” decomposition Q = L^DL (a variant of 
the classical Cholesky decomposition), where L is a lower triangular matrix and D is a 


diagonal matrix. Finally, one obtains r(x) = (L mc;(x))''' D (L md(x)) = Such 

a decomposition is called a sums of squares (SOS) certificate. 


Example 2.2. Let us define r(x) := \ + x\ — 2x\x2 + x^- With m 2 (x) = 

{\,Xi,X 2 ,x\,XiX 2 ,x‘^), one solves the linear matrix feasibility problem r’(x) = 

m 2 (x)''' Qm 2 (x). One can show that the solution writes Q = L^DL for a 6 x 6 matrix 
L and a diagonal matrix D with entries (^, 0,0,1,0, 0), yielding the SOS decomposition: 
r(x) = (i)^ + {xl — X 2 )^. This is enough to prove that p is nonnegative. 


2.2.2. Dense SDP relaxations for polynomial optimization. In order to solve our goal problem 
(Problem 0 )- we are trying to solve Problem ([^, recast as Problem Q. We first explain 
how to obtain tractable approximations of this di fficult problem . Define go := 1. The hier¬ 
archy of SDP relaxations developed by Lasserre [Lasserre 2001 provides lower bounds of 
r*, through solving the optimization problems (Pdp- 


(Pd 


(Pd) : 


sup /i , 

(Tj ,/X 

s.t. r(x) - p = Y!}=o 0'i(x)ffi(x), Vx , 

/r G K, (Tj G E[x], j = 0,..., A:, 
deg(cri5j) < 2d, j = Q,...,k. 


One can solve (Pd) with SDP optimization to find a tuple (/r, ag,..., dk) which enables a 
proof that 51 (x) > 0 A • • • A 5 fc(x) > 0 r(x) — p > 0. 

The next theorem is a consequence of the assumption mentioned in Remark |2.1[ 

Theorem 2.3 (Lasserre [Lasserre 2001| ). Let p^ be the optimal value of the SDP 
relaxation (Pd). Then, the sequence of optimal values (pJ)dGN is nondecreasing and con¬ 
verges to r*. 


The number of SDP variables (i.e. the number of variables of the semidefinite relax¬ 
ation (Pd)) grows polynomially with the integer d, called the relaxation order. Indeed, 
at a fixed number of variables n, the relaxation (Pd) involves 0((2d)"') SDP variables and 
(fc -I- 1) linear matrix inequalities (LMIs) of size 0(d"). When d increases, then more accu¬ 
rate lower bounds of r* can be obtained, at an increasing computational cost. At a fixed 
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d, the relaxation (Pd) involves SDP variables and (d+ 1) linear matrix inequalities 

(LMIs) of size 0{n‘^). 

Example 2.4. Consider the polynomial / mentioned in Section^ fijx) := X 2 X 5 + x^xq — 
X 2 X 3 — x^Xq + Xi{—Xi + X 2 + X 3 — X 4 ^ + X 5 + Xq) and the set K :=^, 6.36]®. The set K can 
be equivalently rewritten as: 


K := {x G 


5i(x) > 0,...,57(x) > 0}, 


with 5 i(x) := (6.36 — Xi){xi — 4) for each i = 1,... ,6 and 57 (x) := 243 — 


Here the 
is 


constant M = 243 is chosen so that M > 6 x 6.36^ and the assumption in Remark 2.1 
fulfilled. The number of initial variables of the optimization problem r* := infxgE" { /(x ) : 
X G K } is n = 6 and the number of optimization constraints is fc = 7. For d = 1, the 
dense SDP relaxation (Pi) involves = 28 variables and provides a lower 

bound Pi = 20.755 for r*. The dense SDP relaxation (P 2 ) involves (^ 4 "^) = 210 variables 
and provides a tighter lower bound of P 2 = 20.8608 for r*. 

2.2.3. Exploiting sparsity. Here we recall how to exploit the structured sparsity of the prob¬ 
lem to replace one SDP problem (Pd) by an SDP problem (Sd) of size where k is 

the average size of the maximal cliques of the correlatio n sparsity pattern (csp) of the poly¬ 
nomial variables (see Waki et al. 2006| Lasserre 2006 for more details). We now present 
these notions as well as the formulation of sparse SDP relaxations (Sd). 

We denote by N" the set of n-tuple of nonnegative integers. The support of a polynomial 
r(x) := defined as supp(r) := {a G N" : ^ 0}. For instance the 

support of r(x) := ^ + xf - 2x1x1+ x^ is supp(p) = { (0, 0), (4, 0), (2,2), (0,4) }. 

Let Fj be the index set of variables which are involved in the polynomial gj, for each 
j = 1,..., /c. The correlative sparsity is represented by the nxn correlation sparsity pattern 
matrix (csp matrix) R defined by: 


R(i, j) := 


1 ifi = j , 

1 if da G supp(/) such that ai,aj > 1, 
1 if 31 G {1, ..., A:} such that i,j G F), 

0 otherwise . 


We define the undirected csp graph G{N, E) with fV = {1,..., n} and E = {{i,j} '■ hj G 
N, i < j, R(i, j) = 1}. Then, let Ci,..., Cm C N denote the maximal cliques of G{N, E) 
and define rij := #61,, for each j = 1 ,..., to. 

Remark 2.5. Assuming that the set K is as in Remark |2.1[ one replaces the constraint 
M — X]r=i ^ 0 by the to redundant additional constraints: 

gk+j ■■= ^ > 0 , j = 1,..., TO , (9) 

i&Ci 

set k' = k + m, define the compact semialgebraic set: 

K' := { X G M” : 51 (x) > 0,..., gfe- (x) > 0 } , 
and modify Problem (|^ into the following optimization problem: 

r* := inf { r(x) : x G K' } . (10) 

xGR" 

For each j = 1,..., to, we note R 2 d[x, Cj] the set of polynomials of K 2 d[x] which involve the 
variables {xi)i^Cj- We denote S[x, Q] := S[x] P| R 2 d[x, Qj. Similarly, we define E[x, F],], 
for each j = 1,..., fc'. The following program is the sparse variant of the SDP program 
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(Pd): 


(Sd) : 


sup /X , 

IJ,,CTj 


s.t. r(x) - fi = , Vx , 

Ai G K , (To G J2]Li S[x, Cj], 

(Tj G S[x, j = 
deg((Tj 5 j) < 2 d, j = 0 ,...,fc', 


where (Tq G if and only if there exist cr^ G S[x, Ci],..., cr"* G S[x, Cm] such 

that cro(x) = J2jLi for all x G M". 

The number of SDP variables of the relaxation (S^) is J2jLi fixed d, it yields 

an SDP problem with 0{k‘^‘^) variables, where k := ^ SjLi i® ffi® average size of the 
cliques Ci,..., Cm- Moreover, the cliques Ci,..., Cm satisfy the running intersection prop¬ 
erty: 


Definition 2.6 {RIP). Let m G No and /i,..., /m be subsets of {1,..., n}. We say that 
/i,..., /m satisfy the running intersection property (RIP) when for all * = 1 ,..., m, there 
exists an integer I < i such that fi n {^j<ilj) C fi. 


This RIP property together with the assumption mentioned in Remark |2 . 5| allow us to state 
the sparse variant of Theorem |2.3| 


Theorem 2.7 (Lasserre [Lasserre 2006 Theorem 3.6]). Let rj he the optimal 


value of the sparse SDP relaxation (S^) 
converges to r* 


Then the sequence {r^)d£f^ is nondecreasing and 


The interested reader can find more details in Waki et al. 2006 about additional ways 


to exploit sparsity in order to derive analogous sparse SDP relaxations. We illustrate the 
benefits of the SDP relaxations (S^j) with the following example: 


Example 2.8. Consider the polynomial / mentioned in Section]^ /(x) := X 2 Xr, + x^XQ — 
X 2 X 3 — x^xe + xi(—xi + X 2 + X 3 — Xi + X 3 + xq). Here, n = 6 ,d = 2, N = {1,..., 6 }. The 
6 x 6 correlative sparsity matrix R is: 


/I 1 1 1 1 1\ 
1110 10 
1110 0 1 
10 0 10 0 
110 0 11 
yi 0 1 0 1 1/ 


The csp graph G associated to R is depicted in Figure [l] The maximal cliques of G are 
Cl := {1,4}, C 2 := {1,2,3}, C 3 := {1,2,5}, C 4 := {1,5,^ and C 5 := {1,3,6}. For d = 2 , 
the dense SDP relaxation (P 2 ) involves (^ 4 "^) = 210 variables against -|- 4 (^ 4 ‘*) = 155 
for the sparse variant (S 2 ). The dense SDP relaxation (P 3 ) involves 924 variables against 
364 for the sparse variant (S 3 ). This difference becomes significant while considering that the 
time complexity of semidefinite progra mming is polynomial w.r.t. the number of variables 
with an exponent greater than 3 (see |Ben-Tal and Nemirovski 2001 Chapter 4] for more 
details). 
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Fig. 1. Correlative sparsity pattern graph for the variables of / from Example |2.8| 


2.3. Computer Proofs for Polynomial Optimization 

Here, we briefly recall some existing features of the COQ proof assistant to handle formal 
polynomial optimization, when using SDP relaxations. The advantage of such relaxations 
is that they provide SOS certificates, which can be formally checked a posteriori. For more 


details on COQ, we recommend the documentation available in Bertot and Casteran 2004 
Given a polynomial r and a set of constraints K, one can obtain a lower bound on r by 
solving any instance of Problem (Pd). Then, one can verify formally the correctness of the 
lower bound rj, using the SOS certificate output cto, ..., <Jk- Indeed it is enough to prove 
the polynomial equality r(x) — rj = inside COQ. Such equalities can be 

efficiently proved using COQ’s ring tactic [Gregoire and Mahboubi 2005 via the mechanism 
of computational reflection [Boutin 1997] . Any polynomial o f type pexprC (see Section m 
can be normalized to a unique polynomial of type polC (see [Gregoire and Mahboubi 2005] 
for more details on the constructors of this type). For the sake of clarity, let us consider 
the unconstrained case, i.e. K = K". One encodes an SOS certificate cro(x) = 9? 

the sequence of polynomials [qi;..qm], each qi being of type polC . To prove the equality 
r = CTo, our version of the ring tactic normalizes both r and the sequence [gi;...; q^] and 
compares the two normalization results. This mechanism is il lustrated in Figure [^ with the 
polynomial r(x) := \ + x\ — 2x\x2 + X 2 (see Example 2.21 being encoded by r and the 
polynomials 1/2 and x\ — X 2 being encoded respectively by qi and q 2 . 


[ I + ^1 - 23:1^2 + xj 


interpretation 


reflexive tactic 




normalization 


reification 


[qcqs] 


Fig. 2. An illustration of computational reflection. 


In the general case, this computational step is done through a checker_sos procedure 
which returns a Boolean value. If this value is true, one applies a correctness lemma, whose 
conclusion yields the nonnegativity of r — over K. In practice, the SDP solvers are 
implemented in floating-point arithmetic, thus the above equality between t — and the 
SOS certificate does not hold. However, following Remark [2.1[ each variable lies in a closed 
interval, thus one can bound the remainder polynomial e(x) := r(x) — i^)9j (^) 

using basic interval arithmetic, so that the lower bound e* of e yields the valid inequality: 
Vx G K,r(x) > -F e*. For more explanation, we refer the interested reader to the formal 
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Input: input variables x, input constraints X, nonlinear expression /, rounded expression 
/, error variables e, error constraints E, relaxation order d 
Output: interval enclosure Id of the error f — f over K X x E 
1: Define the absolute error r(x, e) := /(x,e) — /(x) 

2 : Compute ;(x,e) := r(x, 0) + 0)ej 

3: Define h := r — I 

4: Compute bounds for h: := ia bound(/i, K) 

5: Compute bounds for 1: := sdp_bound(Z, K, d) 

6: return Id := lli + 

Fig. 3. bound: our algorithm to compute roundoff errors bounds of nonlinear programs. 


framework Magron et al. 2015b Section 2.3]. Note that this formal verification remains 
valid when considering the sparse variant (S^). 


3. GUARANTEED ROUNDOFF ERROR BOUNDS USING SDP RELAXATIONS 

In this section, we present our new algorithm, relying on sparse SDP relaxations, t o bo und 
roundoff errors of nonlinear programs. After stating our general alg orith m (Section |3.1[ ), we 
detail how this procedure can handle polynomial programs (Section |3.2[ ) . Extensions to the 
non-polynomial case, including conditional statements, are presented in Section |3.3| 

3.1. The General Optimization Framework 

Here we consider a given program that implements a nonlinear transcendental expression 
/ with input variables x satisfying a set of constraints X. We assume that X is included 
in a box (i.e. a product of closed intervals) [a, b] := [ai,5i] x • • • x [a„,5„] and that X is 
encoded as follows: 


X:={xG 


5 i(x) > 0 , ...,gfe(x) > 0 }, 


for polynomial functions 51 ,..., gu- Then, we de note by /(x, e) the rounded expression of / 
after applying the round procedure (see Section [2.1[ ), introducing additional error variables 
e. ^ 

The algorithm bound, depicted in Figure ^ takes as input x, X, /, /, e as well as the 
set E of bound constraints over e. Here we assume that our program implemen ting / does 
not involve conditional statements (this case will be discussed later in Section 3.31. For a 
given machine e, one has E := [—e, ej™, with m being the number of error variables. This 
algorithm actually relies on the sparse SDP optimization procedure (Sd) (see Section 2.2 
for more details), thus bound also takes as input a relaxation order d G N. The algorithm 
provides as output an interval enclosure Id of the error /(x,e) — /(x) over K. From this 
interval Id '■= [fd, fd], one can compute fd ■= max{—/d, fd}, which is a sound upper bound 
of the maximal absolute error r* := max(x.e)eK I /(x;®) ~ /(^) I- 

After defining the absolute roundoff error r := f — f (Line [^, one decomposes r as the 
sum of an expression I which is affine w.r.t. the error variable e and a remainder h. One way 
to obtain I is to compute the vector of partial derivatives of r w.r.t. e evaluated at (x, 0 ) and 
finally to take the inner product of this vector and e (Line [^. Then, the idea is to compute 
a precise bound of I and a coarse bound of h. The underlying reason is that h involves 
error term products of degree greater than 2 (e.g. 6162 ), yielding an interval enclosure I^ 
of a priori much smaller width, compared to the interval enclosure P of 1. One obtains 

using the procedure ia_bound implementing basic interval arithmetic (Line to bound 

the remainder of the multivariate Taylor expa nsion of r w.r.t. e, expr essed as a combination 
of the second-error derivatives (similar as in Solovyev et al. 2015|). The main algorithm 
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presented in Figure]^ is very similar to the algorithm of FPTaylor Solovyev et al. 2015 


except that SDP based techniques are used instead of the global optimization procedure 


from Solovyev et al. 2015 

. Note that overflow and denormal are neglected here but one 

could handle them, as in 

Solovyev et al. 2015 , by adding additional error variables and 


discarding the related terms using naive interval arithmetic. 


3.2. Polynomial Programs 

We first describe our sdp_bound optimization algorithm when implementing polynomial 

programs. In this case, sdp_bound calls an auxiliary procedure sdp_poly. The bound of I 

is provided through solving two sparse SDP instances of Problem (S^^), at relaxation order 

d. We now give more explanation about the sdp_poly procedure. 

We can map each input variable Xi to the integer i, for all i = 1,..., n, as well as each 
error variable ej to n + j, for all j = 1,..., m. Then, define the sets Ci := {1,... ,n,n + 
1},..., Cm '■= {1, ■ ■ ■ ,n, n + m}. Here, we take advantage of the correlation sparsity pattern 
of I by using m distinct sets of cardinality n +1 rather than a single one of cardinality n + m, 
i.e. the total number of variables. After writing Z(x, e) = r(x, 0) + 

noticing that r(x,0) = /(x, 0) — /(x) = 0, one can scale the optimization problems by 
writing 

771 771 

^(x,e) = ^Sj(x)ej = e^Sj(x)^ , (11) 

3=1 3=1 


with Sj(x) := (x, 0), for all j = 1,..., m. Replacing e by e/e leads to computing an 

interval enclosure of l/e over K' := X x [—1,1]™. Recall that from Remark 2.1 there exists 
an integer M > 0 such that Af —xf > 0, as t he input variables satisfy box constraints. 
Moreover, to fulfil the assumption of Remark |2.5[ one encodes K' as follows: 


K' := { (x,e) e : gi(x) > 0,... ,g/c(x) > 0 , 

5 fc+i(x,ei) > 0 ,...,gk+m{^,em) > 0 }, 


with ej) := M + 1 — ~ J = • ■ • > The index set of variables 

involved in gj is Fj := iV = {1,..., n} for all j = 1,..., k. The index set of variables involved 
in gk+j is Fk+j := Cj for all j = 1 ,..., m. 

Then, one can compute a lower bound of the minimum of r(x,e) := /(x,e)/e = 
Si(x)ej over K' by solving the following optimization problem: 


:= sup g, , 

s.t. k - g = ao + ^393 : 

g€R, 0-0 G E^iT:[(x,e),C'j], 
aj G E[(x,e),Fj], j = 1,.. .,fc + m, 
deg(o'j 5 j) < 2d, j = 1,..., fc + m . 


A feasible solution of Problem (12) ensures the existence of G E[(x, ei)], 
S[(x, Bm)] such that tJo = ’ allowing the following reformulation: 


( 12 ) 


cr™ G 


:= sup g , 

^y(7j 

s.t. I -g = J2 j=i <^393 - 

g GM., Uj G S[x] , j = 1,..., m, 

(jii G S[(x,ej)] ,deg(cr-^) < 2d, j = 1,... ,m, 
deg(crj 5 j) < 2 d, j = l,...,k + m. 


(13) 
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An upper bound l'^ can be obtained by replacing sup with inf and Z' — ^ by ^ — Z' in 
Problem (13l. Our optimization procedure sdp_poly computes the lower bound Z^ as well 


as an upper bound Z^ of Z' over K' then returns the interval := [e Z^, e Z^], which is a sound 
enclosure of the values of Z over K. 

We emphasize two advantages of the decomposition r := I + h and more precisely of 
the linear dependency of I w.r.t. e: scalability and robustness to SDP numerical issues. 
First, no computation is required to determine the correlation sparsity pattern of Z, by 
comparison to the general case. Thu s, it becomes much easier to handle the optimization 
of Z with the sparse SDP Problem (131 rather than with the corresponding instance of 
the dense relaxation (Pd). While the latter involves SDP variables, the former 

involves only m variables, ensuring the scalability of our framework. In addition, 

the linear dependency of I w.r.t. e allows us to scale the error variables and optimize over 
a set of variables lying in K' := X x [—1,1]™. It ensures that the range of input variables 
does not significantly differ from the range of error variables. This conditio n is mandatory 
while considerin g SDP relaxations because most SDP solvers (e.g. Mosek Andersen and 
Andersen 2000 ) are implemented using double precision floating-point. It is impossible to 
optimize Z over K (rather than I' over K') when the maximal value e of error variables is 
less than 2“^^, due to the fact that SDP solvers would treat each error variable term as 0, 
and consequently Z as the zero polynomial. Thus, this decomposition insures our framework 
against numerical issues related to finite-precision implementation of SDP solvers. 

Let us define the interval enclosure P := [Z,Z], with Z := inf(x,e)GK ^(>^^1 e) and Z := 
sup(x e)GK ®)' The next lemma states that one can approximate P as closely as desired 
using the sdp poly procedure. 

Lemma 3.1 (Convergence of the sdp poly procedure). Let P^ he the interval 
enclosure returned by the procedure sdp_poly(Z, K, cZ). The sequence {P^d&i converges to 


Proof. It is sufficient to show the similar convergence result for Z' = Z/e, as it implies 
the convergenc e for I by a scaling argument. The sets Ci,..., Cm satisfy the RIP property 
(see De finiti on 2.6 1 . Moreo ver, the encoding of K' satisfies the assumption mentioned in 

' implies that the sequence of lower bounds (Z(;)deN converges 


Remark 


2.5 


Thus, Theorem 


2.7 


to Z( := inf(-x_e)GK' Z'(x,e). Similarly, the sequence of upper bounds converge to Z', yielding 
the desired result. □ 


Lemma 3.1 guarantees asymptotic convergence to the exact enclosure of I when the re¬ 
laxation order d tends to infinity. However, it is more reasonable in practice to keep 
this order as small as possible to obtain tractable SDP relaxations. Hence, we generi- 
cally solve each instance of Problem ( |13[ ) at the minimal relaxation order, that is do := 
max{ [deg Z/2]), maxi<j <;,+„{ [deg(gj)7^)}}. 


3.3. Non-polynomial and Conditional Programs 

Other classes of programs do not only involve polynomials but also semialgebraic and tran¬ 
scendental functions as well as conditional statements. Such programs are of particular 
interest as they often occur in real-world applications such as biology modeling, space con¬ 
trol or global optimization. We present how the general optimization procedure sdp_bound 

can be extended to these nonlinear programs. 


3.3.1. Semialgebraic programs. Here we assume that the function Z is semialgebraic, that is it 


involves non-polyn omial components such as divisions or square roots. Following Lasserre 
and Putinar 2010 , we explain how to transform the optimization problem inf(x,e)GK Z(x, e) 
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Input: input variables y, input constraints K, semialgebraic expression / 
Output: variables Ypoiy, constraints Kpoiy, polynomial expression /poly 
1: I := ia bound(/, K) 

2: if / = Pol (p) then ypoiy := y, Kpoiy := K, /poiy := p 


3 

else if / 

= Div 

(g,h) then 


4 


ffpoly : = 

lift(y,K,g) 


5 

y/i,K,,, 

^poly .— 

lift(y,K,/i) 


6 

ypoiy := 

(yg^Yh 

2:) /poly 

:= X 

7 

Kpoiy 

= {ypoiy 

G Kg X K/, X I 

■ ^^poly 

8 

else if / 

= Sqrt 

(g) then 


9 

yg, Kg, 

ffpoly : = 

lift(y,K,5) 


10 

: ypoiy : = 

= (yg, 2 :) 

/poly : = 

X 

11 

: Kpoiy : 

II 

G Kg X I: x^ = 

- 5 poly} 

12 





13 

: end 




14 

: return ypoiy 5 ^poly;/poly 



Fig. 4. lift: a recursive procedure to reduce semialgebraic problems to polynomial problems. 


into a polynomial optimization problem, then use the sparse SDP program (13l. One way 
to perform this reformulation consists of introducing lifting variables to represent non¬ 
polynomial operations. We first illustrate the extension to semialgebraic programs with an 
example. 


Example 3.2. Let us consider the program implementing the rational function / : [0, 1] —>■ 
K defined by f{xi) := Applying the rounding procedure (with machine e) yields 

f{xi,e) := and the decomposition r(a;i,e) := f{xi,e) - f{xi) = l{xi,e) + 

h{xi,e) = si{xi)ei + 52 ( 2 : 1)62 -F h{xi,e). One has Si(a:i) = (a:i, 0) = and 

52 ( 2 : 1 ) = -si(a:i). 

Let K := [0, 1] X [—e, e]^. One introduces a lifting variable X 2 '■= to handle the 
division operator and encode the equality constraint p(x) := 0 : 2(1 + 2 : 1 ) — a;i = 0 with the 
two inequality constraints p(x) > 0 and — p(x) > 0. To ensure the compactness assumption, 
one bounds X 2 within / := [0, 1 / 2 ], using basic interval arithmetic. 

Let Kpoiy := {(x,e) e [0, 1] x / x [e,e]^ : p(x) > 0, —p(x) > 0}. Then the rational 
optimization problem involving I is equivalent to inf(x,e)GKpoiy 2 : 2 (—ei + 62 ), a polynomial 
optimization problem that we can handle with the sdp_poly procedure, described in Sec¬ 

tion |221 


In the semialgebraic case, sdp_bound calls an auxiliary procedure sdp_sa. Given input 

variables y (x,e), input constraints K := X x E and a semialgebraic function I, sdp_sa 

first applies a recursive procedure lift which returns variables ypoiy, constraints Kpoiy and 
a polynomial /poly such that the interval enclosure P of l{y) over K is equal to the inter¬ 
val enclosure of the polynomial ^poiy(ypoiy) over Kpoiy. Calling sdp_sa yields the interval 
enclosure := sdp_poly(Zpoi;j^Kpoiy, d). We detail the lifting procedure lift in Figured 
for the constructors Pol (Line |^, Div (Line and Sqrt (Line |^. The interval / obtainea 
through the ia_bound procedure (Line (ly) fflows us to constrain the additional variable x 


to ensure the assumption of Remark [2.5| For the sake of consistency, we omit the other cases 
(Neg, Add, Mul and Sub) where the procedure is straightforward. Fo r a similar pro cedure in 
the context of global optimization, we refer the interested reader to Magron 2013 Chapter 
2]. The set of variables ypoiy can be decomposed as (xpoiy,e), where Xpoiy gathers input 
variables with lifting variables and has a cardinality equal to ripoiy. Then, one easily shows 
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that the sets ripoiy, ei},. • ■ ,{1, ■ • •, «poiy, em} satisfy the RIP, thus ensuring to solve 

efficiently the corresponding instances of Problem (Hi. 


3.3.2. Transcendental programs. The above lifting procedure allows an exact representation 
of the graph of a semialgebraic function with polynomials involving additional (lifting) vari¬ 
ables. We consider a procedure to approximate transcendental functions with semialgebraic 
functions. Here we assume that the function I is transcendental, i.e. involves u nivariate 
non-semialge braic components such as exp or sin. We use the method presented in [Magron 


et al. 2015a , based on maxplus approximation of semiconvex transcendental functions by 


quadratic functions. This idea comes from optimal control McEneaney 2006 and was de¬ 


veloped further to represent the value function by a “maxplus linear combination”, which 
is a supremum of quadratic polynomials at given points Xi. Given a set of points {xi), we 
approximate from above and from below every transcendental function /r by infima and 
suprema of finitely many quadratic polynomials (/“) and Hence, we reduce the prob¬ 

lem to semialgebraic optimization problems. We can interpret this method in a geometrical 
way by thinking of it in terms of “quadratic cuts”, since quadratic inequalities are added to 
approximate the graph of a transcendental function. 

For each univariate transcendental function /r in our dictionary set T>, one assumes that 
/r is twice differentiable, so that the univariate functio n fl := /r -I- ^ | ■ P i s convex on I for 
large enough 7 > 0 (for more details, see the reference McEneaney 2006 ). It follows that 
there exists a constant 7 < sup^jg/ such that for all Xi £ 1: 


Vxe/, fuix) > f^^{x), 

- 7 2 (14) 

with f~ := -T(^x- XiY + f^{xi){x - Xi) + fR{xi) , 

implying that for all x € I, fR{x) > max^^^g//“(a;). Similarly, one obtains an upper- 
approximation mina;^g 7 /+(a:). Figure pi provides such approximations for the function 
/mix) := log(l -|-exp(a:)) on the intervalT := [— 8 , 8 ]. 



Fig. 5. Semialgebraic Approximations for x 1 —> log(l + exp(x)): max{fg {x),fg^ (x)} < log(l + exp(a;)) < 
min{/o+(a:),/+(a;)}. 


For transcendental programs, our procedure sdp_bound calls the auxiliary procedure 

sdp_transc. Given input variables (x,e), constraints K and a transcendental function I, 

sdp_transc first computes a semialgebraic lower (resp. upper) approximation l~ (resp. fy) 

of I over K. For more d etails in the context of global optimization, we refer the reader 
to [Magron et al. 2015a . Then, calling the procedure sdp_sa allows us to get interval 
enclosures of L~ as well as fy. We illustrate the procedure to handle transcendental programs 
with an example. 
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Input: input variables x, input constraints X, nonlinear expression /, rounded expression 
/, error variables e, error constraints E, relaxation order d 
Output: interval enclosure Id of the error f — f over K X x E 
1: if / = IfThenElse (p,g^h) then 
2 : := bound(x,X,p,p,e, E, d) = [pd,Pd] 

3 : Xi := {x G X : 0 < p(x) < -pd} 

4 : X2 := {x G X : -p^ < p(x) < 0 } 

5: X3 := {x G X : 0 <p(x)} 

6: X4 := {x G X : p(x) < 0 } 

7 : /i := bound_nlprog(x, Xi, g, /i, e, E, d) 

8: /j := bound_nlprog(x, X2,/i, 5, e, E, d) 

9 : ij := bound_nlprog(x, X3, g, 5, e, E, d) 

10: := bound_nlprog(x, X4, h, fi, e, E, d) 

11: return Id := I^UljUl^Ulj 

12: else return Id := bound(x, X, /,/, e, E, d) 

13 : end 

Fig. 6. bound_nlprog: our algorithm to compute roundoff error bounds of programs with conditional state¬ 
ments. 


Example 3 . 3 . Let us consider the program implementing the transcendental function 
/ : [—8,8] —> K defined by f{xi) := log(l + exp(a:i)). Applying the rounding procedure 
yields f{xi,e) := log[(l + exp(a;i)(l + ei)) (1 + e2)](l + 63). Here, |e2| is bounded by the 
machine e while jeij (resp. [631) is bounded with an adjusted absolute error ei := e(exp) 
(resp. 63 := e(log)). Let K := [-8,8] x [-ei,ei] x^[-e,e] x [-e3,e3]- 

One obtains the decomposition r{xi,e) := /(a;i,e) — f{xi) = l{xi,e) + h{xi,e) = 
si{xi)ei + S2{xi)e2 + S3(a;i)e3 + h{xi,e), with si(a:i) = S2{xi) = 1 and 

53(0:1) = log(l + exp(a;i)) = f{xi). Figure @ provides a lower approximation Sg := 
max{/^,/g"} of S3 as well as an upper approximation s]]" := min{/^,/^}. One can get 
similar approximations sj” and sj*" for si. One first obtains (coarse) interval enclosures 

I2 = ia bound(si,K) and J3 = ia bound(s3,K) and one introduces extra variables 

0:2 G I2 and 0:3 G /s to represent si and S3 respectively. Then, the interval enclosure of 
I over K is equal to the interval enclosure of 4a(x:,e) := 0:261 + 62 + 0:363 over the set 
Ksa := {(a:i,e) G K , (0:2,0:3) G h x h ,Si{xi) < X2 < sf{xi) ,S;]'(o;i) < 0:3 < s;]'(a;i)}. 


3.3.3. Programs with conditionals. Finally, we explain how to extend our bounding proce¬ 
dure to nonlinear programs involving conditionals through the recursive algorithm given 
in Fi gure |6l The overall procedur e is very similar to the one implemented within the Rosa 
tool Darulova and Kuncak 2014| Section 7, Figu re 6 ], The bound nlprog algorithm relies 
on the bound procedure (see Figure in Section |3.1[ ) to compute roundoff error bounds of 
programs implementing transcendental functions (Line From Line rm 


algorithm handles the case when the program implements a function / d^ned as follows: 


a to Line 
ned as follo' 


the 


jl-x) ■= ifp(x)>0 , 

\ft.(x) otherwise. 


The first branch output is g while the second one is h. More sophisticated conditionals, such 
as “_pi(a:) > 0 or/and P 2 {x) > 0”, are not handled at the moment but one could easily extend 
the current framework to do so. A preliminary step consists of computing the roundoff error 
enclosure I^ := [i^,Pd] (Line 0 for the program implementing the polynomial p. Then the 


ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY. 











A:19 


procedure computes bounds related to the discontinuity error of the branch, that is the 
maximal value between the four following errors: 

— (Line the error obtained while computing the rounded result h of the second 
branch instead of computing the exact result g of the first one, occurring for the set 
of variables (x,e) such that p(x,e) < 0 < p(x). For scalability and numerical is¬ 
sues, we consider an over-approximation Xi (Line of this set, where the variables 
X satisfy the relaxed constraints 0 < p(x) < —pd- Note that in this case, one has 


Z(x,e) = r(x, 0 ) + Er=i 


m ^^(Xje) , 


^j=i with r(x,0) = h(x) - g{x) ^ 0. In general, 

we expect the magnitude of the partial derivative sum to be very small compared to the 
one of r(x, 0 ). 

(Line ^ the error obtained while computing the rounded result g of the first branch 
insteacl of computing the exact result h of the second one, occurring for the set of 
variables (x,e) such that p(x) < 0 < p(x,e). We also consider an over-approximation 
X 2 (Line 0, where the variables x satisfy the relaxed constraints —pd < _p(x) < 0. 
the roundoff error corresponding to the program implementation of g. 
the roundoff error corresponding to the program implementation of h. 



3.3.4. Simplification of error terms. In addition, our algorithm bound_nlprog integrates sev¬ 
eral features to reduce the number of error variables. First, it memorizes all sub-expressions 
of the nonlinear expression tree to perform common sub-expressions elimination. We can 
also simplify error term products, thanks to the following lemma. 


Lemma 3.4 (Higham |Higham 2002 Lemma 3.3]). 
and assume that for a given integer k, one has e < ^ 
ei,..., Cfc G [—e, e], there exists 9k such that nti(l + e,) 


Let e be the machine precision 
and 7 fc := Then, for all 

= l + 9k and \ 9k \< 7fc. 


Lemma 


3.4 


implies that for any k such that e < ^, one has 9k < {k+ l)e. Our algorithm has 
an option to automatically derive safe over-approximations of the absolute roundoff error 
while introducing only one variable ei (bounded by (fc -|- l)e) instead of k error variables 
ei,..., Cfc (bounded by e). The cost of solving the corresponding optimization problem can 
be significantly reduced but it yields coarser error bounds. 


4. EXPERIMENTAL EVALUATION 

Now, we present experimental results obtained by applying our general bound_nlprog 

algorithm (see Section]^ Figure]^ to various examples coming from physics, biology, space 

control and optimization. The bound_nlprog algorithm is implemented in an open-source 

tool called Real2Float, built in top of the NLCertify nonlinear verification package, relying 
on OCaml (Version 4.02.1), COQ (Version 8.4pl5) and interfaced with the SDP solver Sdpa 
(Version 7.3.9). The SDP solver output numerical SOS certificates, which are converted 
into rational SOS using the Zarith OCaml library (Version 1.2), implementing arithmetic 
operations over arbitrary-precision integers. For more details about the installation and 
usage of Real2Float, we refer to the dedicated web-pag^and the setup instructions]^ All 
examples are displayed in Appendix A as the correspo nding Real2Float input text files and 
satisfy our nonlinear program semantics (see Section HI- All results have been obtained 
on an Intel Core i7-5600U CPU (2.60GHz). Execution timings have been computed by 
averaging over five runs. 


" http://nl-certify.forge.ocamlcore.org/real2float.html 
^see the README. md file in the top level directory 


ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY. 











A:20 


4.1. Benchmark Presentation 

For each example, we compared th e qu ality of the roundoff error bounds (Table |n| and cor¬ 
responding executio n times (Table |III|) wh ile running our tool Real2Float, FPTa ylor (ver¬ 
sion from May 201 6 Solovyev et ai. 2(jl5 ), Rosa (version from May 2014 used in |Darulova 
and Kuncak 2014 ), Gappa (version 1.2.0 Daumas and Melquiond 2010 ) and Fluctuat 
(version 3.1370 Delmas et al. 2009 ). 


To ensure fair comparison, our initial choice was to focus on tools providing certificates 
and using the same rounding model (FPTaylor or Rosa which relies on an SMT solver 
theoretically able to output satisfiability certificates). However, for the sake of complete- 


ness, we have also compared Real2Float with Gappa Daumas and Melquiond 2010 and 


Fluctuat Delmas et al. 2009 . For all tools, we use default parameters, we use the default 


number of subdivisions in Fluctuat and for Gappa, we provide the simplest user-provided 
hints we could think of. 

A head-to-head comparison is not straightforward here due to differences in the ap¬ 
proaches: Gappa uses a n im proved rounding model based on a piecewise constant absolute 
error bound (see Section [T^ for more details), and Fluctuat does not produce output cer¬ 
tificates. We also performed further experiments while turning on the improved rounding 
model of FPTaylor (which is the same as in Gappa and Fluctuat). 

A given program implements a nonlinear function /(x), involving variables x lying in a 
set X contained in a box [a,b]. Applying our rounding model on / yields the nonlinear 
expression /(x,e), involving additional error variables e lying in a set E. 

At a given semidefinite relaxation order d, our tool computes the upper bound fd of 
the absolute roundoff error | / — / | over K := X x E and verifies that it is less than a 
requested number e^eai 2 Fioat • keep the relaxation order d as low as possible to ensure 

tractable SDP programs, it can happen that fd > e^eai 2 Fioaf The Real2Float tool has 
the default option to perform box subdivisions when the number of initial variables and 
maximal polynomial degree are both small. When the option is disabled, the solver does 
not perform subdivisions and outputs the error bound fd- When enabled, we subdivide a 
randomly chosen interval of the box [a, b] in two halves to obtain two sub-sets Xi and X 2 , 
fulfilling X Xi U X 2 , and apply the bound_nlprog algorithm on both sub-sets either 


until we succeed to certify that 


'Real2Float 


is a sound upper bound of the roundoff error or 


until the maximal number of branch and bound iterations is reached. For each benchmark. 


an error bound eReai 2 Fioat is automatically computed while setting e, 


Real2Float 


= 0 . 


The number eReai 2 Fioat is compared with the upper bounds computed by two other tools 
implementing simple rounding models: FPTaylor, which relies on Taylor Symbolic expan- 


sions 


and Kuncak 2014 


Solovyev et al. 2015 , and Rosa, which relies on SMT and affine arithmetic Darulova 


F'or comparison purpose, we also executed each program using random inputs, following 


the approach used in the Rosa paper Darulova and Kuncak 2014 . Specifically, we executed 


each program on 10^ random inputs satisfying the input restrictions. The results from these 
random samples provide lower bounds on the absolute error. For consistency of compar¬ 
ison, the_erro£_bounds_computed with FPTaylor correspond to the procedure FPT. (a) 
(see [Solovy^ et al. 20T5 ) using the same simplified rounding model as the one described 
in Equation |3l also used in Rosa Darulova and Kuncak 2014]. For the sake of further pre¬ 
sentation, we identify with a letter (from a to z and from a to y) ea ch of the 30 nonlinear 
progr ams. The programs a-b and f-s are tak en from the Rosa pap er Darulova and Kuncak 
2014 and were used in the FPTaylor paper Solovyev et al. 2015| as well: 


The first 9 programs implement polynomial functions: a-b come from physics. 


derived from expressions involved in the proof of Kepler Conjecture Hales 2006 


e are 
and 


f-h implement polynomial approximations of the sine and square root functions. The 
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Table II. Comparison results of upper and lower bounds for absolute roundoff errors among tools implementing either simple 
or advanced rounding model. For each model, results of the winning tool are emphasized using bold fonts. 


Benchmark 

id 

Simple rounding 

Real2Float Rosa FPTaylor 

Improved rounding 

FPTaylor GappA FluCTUAT 

lower bound 

Programs involving polynomial functions 

rieidBodvl 

a 

5.33e-13 

5.08e-13 

3.87e-13 

2.95e-13 

2.95e-13 

3.22e-13 

2.28e-13 

rigidBody2 

b 

6.48e-ll 

6.48e-ll 

5.24e-ll 

3.61e-ll 

3.61e-ll 

3.65e-ll 

2.19(^11 

keplerO 

c 

1.18e-13 

1.16e-13 

1.05e-13 

7.47e-14 

1.12(^13 

1.26e-13 

2.23e-14 

keplerl 

d 

4.47e-13 

6.49e-13 

4.49e-13 

2.87e-13 

4.89(^13 

5.57e-13 

7.58e-14 

kepler2 

e 

2.09e-12 

2.89e-12 

2.10e-12 

1.58e-12 

2.45(^12 

2.90e-12 

3.03(^13 

sineTaylor 

f 

6.03e-16 

9.56e-16 

6.75e-16 

4.44e-16 

8.33e^02 

6.86e-16 

2.85e-16 

sineOrderS 

g 

1.19e-15 

l.lle-15 

9.97e-16 

7.95e-16 

7.62e-16 

1.03e-15 

3.34e-16 

sqroot 

h 

1.29e-15 

8.41e-16 

7.13e-16 

5.02e-16 

5.37e^l6 

3.21e-13 

4.45e-16 

himmilbeau 

i 

1.43e-12 

1.43e-12 

1.32e-12 

1 .Ole-12 

1 .Ole-12 

1 .Ole-12 

1.47e-13 

Programs involving semialgebraic functions 

dopplerl 

j 

7.65e-12 

4.92e-13 

1.59e-13 

1.29e-13 

1.82e-13 

1.34e-13 

7.11e-14 

doppler2 

k 

1.57e-ll 

1.29e-12 

2.90e-13 

2.39e-13 

3.23(^13 

2.53e-13 

1.14(^13 

dopplerS 

1 

8.55e-12 

2.03e-13 

8.22e-14 

6.96e-14 

9.29e-14 

7.36e-14 

4.27e-14 

verhulst 

m 

4.67e-16 

6.82e-16 

3.53e-16 

2.50e-16 

3.18e^l6 

4.84e-16 

2.23e-16 

carbonGas 

n 

2.21e-08 

4.64e-08 

1.23e-08 

7.77e-09 

8.85e^09 

1.86e-08 

4.11e-09 

predPrey 

0 

2.52e-16 

2.94e-16 

1.89e-16 

1.60e-16 

1.95e^l6 

2.45e-16 

1.47e-16 

turbine1 

p 

2.45e-ll 

1.25e-13 

2.33e-14 

1.67e-14 

3.88e^l4 

6.09e-14 

1.07e-14 

turbine2 

q 

2.08e-12 

1.76e-13 

3.14e-14 

2.Ole-14 

3.97e^l4 

8.96e-14 

1.43e-14 

turbines 

r 

l.Tlc-ll 

8.50e-14 

1.70e-14 

9.58e-15 

9.96e+00 

4.90e-14 

5.33e-15 

jetEngine 

s 

OoM 

1.62e-08 

1.50e-ll 

1.03e-ll 

1.32e+05 

1.82e-ll 

5.46e-12 


Programs implementing polynomial functions with polynomial preconditions 


floudas2_6 

t 

5.15e-13 

5.87e-13 

7.88e-13 

5.94e-13 

5.98(^13 

7.45e-13 

4.56(^14 

floudas3_3 

u 

5.81e-13 

4.05e-13 

5.76e-13 

4.29e-13 

2.65e-13 

4.32e-13 

1.48e-13 

floudas3_4 

V 

2.78e-15 

2.56e-15 

2.23e-15 

1.78e-15 

1.23e-15 

2.23e-15 

3.80e-16 

floudas4_6 

w 

1.82e-15 

1.33e-15 

1.23e-15 

8.89e-16 

8.89e-16 

1.12e-15 

2.35e-16 

floudas4 7 

X 

1.06e-14 

1.31e-14 

1.80e-14 

1.32e-14 

7.44e-15 

1.71e-14 

7.31e-15 

Programs involving conditional statements 

cavlO 

y 

2.91e+00 

2.91e+00 

- 

- 

- 

1 .02e+02 

2.90e+00 

perin 

z 

2 .01e+00 

2 .01e+00 

- 

- 

- 

4.91e+01 

2.00e+00 

Programs implementing transcendental functions 

logexp 

a 

2.52e-15 

- 

2.07e-15 

1.99e-15 

- 

- 

1.19(^15 

sphere 

p 

1.53e-14 

- 

1.29e-14 

8.21e-15 

- 


5.05e-15 

hartmanS 

Y 

2.99e-13 

- 

1.34e-14 

4.97e-15 

- 

- 

l.lOe-15 

hartmanG 

6 

5.09e-13 

- 

2.55e-14 

8.19e-15 

- 

- 

2.20e-15 


program i is is sued from the gl obal optimization literature and implements the problem 
Himmilbeau in Ali et al. 2005 . 

The 10 programs j-s implement semialgebraic functions: j-1 and p-s come from physics, 
m and o from biology, n from cont rol. All these programs are used to compare FPTaylor 
and Rosa in Solovyev et al. 2015 . 

X come from the global optimization literature and correspond 

We 


— The five programs t 


respectively to Problem 2.6, 3.3, 3.4, 4.6 and 4.7 in iFloudas and Pardalos 1990 

or 


selected them as they typically involve nontrivial polynomial preconditions (i.e. X is 
not a simple box but rather a set defined with conjunction of nonlinear polynomial 
inequalities). 

The two programs y-z involve conditional statements and come from the static analysis 


litera ture. They correspond to the two respective runn ing examples of Ghorbal et al. 
2010| (Fluctuat’s divergence error comp utation) and Alexandre Marechal 2014 . 'Fhe 
first program y is used in the Rosa paper Darulova and Kuncak 2014 for the analysis 
of branches discontinuity error. 

The last four programs a-y invo lve transcendental fu nctions. The two programs a and p 
are used in the FPTaylor paper [Solovyev et al. 2015 and correspond respectively to the 
program logexp (see Example |3.3| and the prograin sphere taken from NASA World 
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Wind Java SDK Nasa 2011 . The 2 programs y and B respectively implement the func¬ 


tions coming from the optimization problems Hartman 3 and Hartman 6 in Ali et al. 


2005 , involving both sums of exponential functions composed with quadratic polynomi¬ 
als. 


Tool comparison settings. The five tools Real2Float, Rosa, FPTaylor, Gappa and Fluc- 
TUAT can handle programs with input variable uncertainties as well as any floating-point 
precision, but for the sake of conciseness, we only considered to compare their performance 
on programs implemented in double precision floating-point (e = 2“®^). By contrast with 
preliminary experiments presented in Section o where we considered floating-point input 
variables, we run each tool by considering all input variables as real variables, thus we apply 
the rounding operation to all of them. For the programs involving transcendental functions, 
we followed the same procedure as in FPTaylor while adjusting the precision e (/h) = 1.5e 
for each special function /r S {sin, cos, log, exp}. Each univariate transcendental function 
is approximated from below (resp. from above) using suprema (resp. infima) of linear or 
quadratic polynomials (see Example |3 .3 1 for the case of program logexp). 


4.2. Comparison Results 

Comparison results for error bound computation are presented in Table [IT] Among the tools 
using the simple rounding model (Real2Float, Rosa and FPTaylor), our Real2Float tool 
computes the tightest upper bounds for 7 (resp. 4) out of 30 benchmarks when comparing 
with Rosa and FPTaylor (resp. all tools). For all programs j-s involving semialgebraic 
functions and the four programs ol-8 involving transcendental functions, FPTaylor computes 
the tightest bounds, when comparing with Real2Float and Rosaj^One current limitation 
of Real2Float is its limited ability to manipulate symbolic expressions, e.g. computing 
rational function derivatives or yielding reduction to the same denominator. In particular, 
the analysis of program s aborted after running out of memory (meaning of the symbol 
OoM). The FPTaylor tool is better suited to handle programs that exhibit such rational 


funct ions and also includes an interface with the Maxima computer algebra system Maxima 
2013] to perform symbolic simplifications. 


We mention that the interested reader can find more detailed experimental comparisons 
between the th ree tools implementi ng improved rounding models (FPTaylor, Fluctuat 
and Gappa) in Solovyev et al. 2015 Section 5.2]. Note that FPTaylor provides the tightest 
bounds for 24 out of 30 benchmarks. The Gappa software provides the tightest bounds for 
8 o ut of 30 benchmarks while being almost always faster than other tools. As emphasized 

Gappa can sometimes compute more precise bounds with more 


Solovyev et al. 2015 


advanced user-provided hints. Note that Fluctuat provides tighter bounds than Gappa 
for most rational functions while being always slower. It would be worth implementing the 
same improved rounding model in Real2Float to perform numerical comparisons. 

To the best of our knowledge, Real2Float is the only academic tool which is able to 
handle the general class of programs involving either transcendental functions or condi¬ 
tional statements. The FPTaylor (resp. Rosa) tool does not currently handle conditionals 
(resp. transcendental functions), as meant by the symbol — in the corresponding column 
entries. However, an interface bridging the FPTaylor and Rosa tools would provide each 
other with the relevant missing features. These error boun d co mparison results together 
with their correspondiM execution timings (given in Table III I are used to plot the data 
points s hown in Figure^ For each benchmark identified by id, let tReai 2 Fioat (in 3rd column 
of Table III I refer to the execution time of Real2Float to obtain the corresponding upper 
bound eReai 2 Fioat (in 3rd column of Table O. 


®The running execution times of Rosa may change with more recent versions 
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Table III. Comparison of execution times (in seconds) for absolute roundoff error bounds among tools im¬ 
plementing either simple or advanced rounding model. For each model, the winner results are emphasized 
using bold fonts. 


Benchmark 

id 

Simple rounding 
Real2Float Rosa FPTaylor 

Improved rounding 
FPTaylor GappA FluCTUAT 

rigidBodyl 

a 

0.58 

0.13 

1.84 

0.41 

0.10 

0.29 

rigidBody2 

b 

0.26 

2.17 

3.01 

0.46 

0.15 

0.52 

keplerO 

c 

0.22 

3.78 

4.93 

6.96 

0.44 

0.70 

keplerl 

d 

17.6 

63.1 

9.33 

4.90 

0.72 

0.37 

kepler2 

e 

16.5 

106 

19.1 

13.8 

1.58 

2.04 

sineTaylor 

f 

1.05 

3.50 

2.91 

0.11 

0.16 

0.55 

sineOrderS 

g 

0.40 

0.48 

1.90 

0.40 

0.06 

0.22 

sqroot 

h 

0.14 

0.77 

2.70 

0.44 

0.19 

0.40 

hiininilbeau 

i 

0.20 

2.51 

3.28 

0.49 

0.09 

1.00 

dopplerl 

j 

6.80 

6.35 

6.13 

1.34 

0.08 

0.77 

doppler2 

k 

6.96 

6.54 

6.88 

1.57 

0.08 

0.79 

dopplerS 

1 

6.84 

6.37 

9.13 

1.50 

0.07 

0.78 

verhulst 

m 

0.51 

1.36 

1.37 

0.40 

0.05 

0.25 

carbonGas 

n 

0.83 

6.59 

3.73 

0.52 

0.15 

2.03 

predPrey 

0 

0.87 

4.12 

1.78 

0.48 

0.04 

5.14 

turbinel 

p 

72.2 

3.09 

4.38 

0.53 

0.21 

5.79 

turbine2 

q 

4.72 

7.75 

3.25 

0.60 

0.12 

4.76 

turbines 

r 

74.5 

4.57 

3.46 

0.61 

0.20 

5.84 

jetEngine 

s 

OoM 

125 

9.79 

3.06 

0.31 

31.2 

floudas2_6 

t 

2.49 

159 

15.9 

14.3 

2.35 

26.8 

floudas3_3 

u 

0.45 

13.9 

5.64 

13.9 

0.76 

6.41 

floudas3_4 

V 

0.09 

0.49 

1.47 

1.27 

0.07 

1.01 

floudas4_6 

w 

0.07 

1.20 

0.91 

0.37 

0.05 

1.69 

floudas4 7 

X 

0.13 

21.8 

1.64 

0.39 

0.07 

0.53 

cavlO 

y 

0.23 

0.59 

- 

- 

- 

1.26 

perin 

Z 

0.49 

2.74 

- 

_ 

- 

1.19 

logexp 

a 

1.05 

- 

1.10 

0.39 

- 

- 

sphere 

p 

0.05 

- 

2.04 

3.69 

- 

- 

hartmanS 

Y 

2.02 

- 

32.5 

27.8 

- 

- 

hartmanB 

6 

119 

- 

364 

259 

- 

- 


Now, let us define the execution times tRosa, tppTayior and the corresponding error bo unds 
CRosa, CFPTayior- Then the x-axis coordinate of the point (resp. [^) displayed in Figure 7(a) 


(resp. |7(b)| ) corresponds to the logarithm of the ratio between the execution time of Ro 
(resp. FPTaylor) and Real2Float, i.e. log — (resp. log ). Similarly, the y-axis 


coordinate of the point (id) (resp. lidl) is log —— (resp. log 

^ N—y V ^ I 1'' ^ eReal2Float ' ^ ^ eReal2Float ' 

The axes of the coordinate system within Figure 0 divide the plane into four quad¬ 
rants: the nonnegative quadrant (-F, -F) contains data points referring to programs for which 
Real2Float computes the tighter bounds in less time, the second one (-F, —) contains points 
referring to programs for which Real2Float is faster but less accurate, the non-positive 
quadrant (—, —) for which Real2Float is slower and computes coarser bounds and the last 
one (—,-F) for which Real2Float is slower but more accurate. 

On the quadrant (-F, —) of Figure [7(b)j one can see that Real2Float computes bounds 
which are less accurate than FPTaylor on semialgebraic and transcendental programs, but 
does so more quickly for most of them. The quadrant (—,—) indicates that Rosa and 
FPTaylor are more precise and efficient than Real2Float on the thr ee pr ograms p-r. The 
presence of 18 plots on the nonnegative quadrant (-F,-F) of Figure [7(^ and Figure [7(b)] 
confirms that Real2Float does not compromise efficiency at the expense of accuracy, in 
particular for programs implementing polynomials with nontrivial polynomial precondi¬ 
tions. 
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For each program implementing polynomials, our tool has an option to provide for¬ 
mal guarantees for the corresponding roundoff error bound eReai 2 Fioat- Using the formal 
mechanism described in Section [231 Real2Float formally checks inside COQ the SOS cer¬ 
tificates generated by the SDP solver for interval enclosures of linear error terms 1. The 
FPTaylor (resp. Gappa) software has a similar option to provide formal scripts which can 
be checked inside the Hol-light (resp. Coq) proof assistant, for programs involving poly¬ 
nomial and transcendental functions. For formal verification, our tool is limited compared 
with FPTaylor and Gappa as it cannot handle non-polynomial programs. 


Table IV. Comparisons of informal and formal execution times to certify roundoff error bounds obtained 
with Real2Float, FPTaylor and Gappa. 


Benchmark 

id 

Informal execution time 
Real2Float FPTaylor GappA 

Formal execution time 
Real2Float FPTaylor GappA 

rigidBodyl 

a 

0.58 

1.84 

0.10 

0.36 

10.2 

2.58 

rigidBody2 

b 

0.26 

3.01 

0.15 

4.81 

32.3 

4.90 

keplerO 

c 

0.22 

4.93 

0.44 

0.29 

45.5 

13.0 

keplerl 

d 

17.6 

9.33 

0.72 

449 

90.5 

28.9 

kepler2 

e 

16.5 

19.1 

1.58 

297 

274 

59.3 

sineTaylor 

f 

1.05 

2.91 

0.16 

7.54 

42.1 

13.6 

sineOrderS 

g 

0.40 

1.90 

0.06 

0.34 

10.4 

6.74 

sqroot 

h 

0.14 

2.70 

0.19 

0.80 

16.8 

12.9 

himmilbeau 

i 

0.20 

3.28 

0.09 

0.89 

32.2 

3.73 


Next, we describe the formal proof results obtained while verifying the bounds for the nine 
polynomial programs a-i. Table [TVl allows us to compare the execution times of Real2Float, 
FPTaylor and Gappa required to analyze the nine programs in both informal (i.e. without 
verification inside GOQ or Hol-light) and formal settings. Comparing Real2Float with 
FPTaylor, the latter performs better than the former to analyze formally the two programs 
d-e but yields coarser bounds. Our formal procedure is less computationally efficient when 
the degree (resp. number of variables) of the SOS polynomials grows, as for these two 
programs. The informal speedup ratio is greater than 3 for the five programs a, c and g-i. 
The table shows that the speedup ratio in the formal setting is higher than in the informal 
setting for these five programs. The explanation could be that GOQ is inherently faster than 
Hol-light at checking computations while delegating expensive ones to a virtual machine 
(being part of Goq’s trusted base). Either SOS or Taylor methods could work with the two 
proof assistants Hol-light and GOQ (with natural changes in execution time) and it would 
be meaningful to compare similar methods in a given proof assistant (SOS techniques in 
Hol-light or Taylor methods in Goq) but both are not yet fully available. The Gappa 
tool performs better when analyzing the two programs d-e, but is less accurate. This is in 
contrast with all other programs where the formal verification with Real2Float is faster 
than with Gappa while providing coarser bounds (except program f where Gappa yields 
pessimistic results) Our formal verification framework is a work in progress as it only allows 
to check correctness of SOS certificates for polynomial programs. The step from formal 
verification of the step translating the inequality to a statement of explicit syntactic form 
“Vx, |rounded(/(x)) — /(x)| < e” is missing and requires more software engineering. Thus, 
the timings presented in Table IV for Real2Float do not include formal computation of the 
symbolic second-order derivatives of r w.r.t. e as well as the cost of bounding them using 
formal interval arithmetic. 

At the moment (and in contrast to our method) the Rosa tool does not formally verify the 
output bound provided by the SMT solve r, but such a feature could be embedded through 
an interface with the smtcoq framework |Armand et al. 201 1|. This latter tool allows the 
proof witness generated by an SMT solver to be formally (arid independently) re-checked 
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inside COQ. The smtcoq framework uses tactics based on computational reflection to enable 
this re-checking to be performed efficiently. 

5. CONCLUSION AND PERSPECTIVES 

Our verification framework allows us to over-approximate roundoff errors occurring while 
executing nonlinear programs implemented with finite precision. The framework relies 
on semidefinite optimization, ensuring certified approximations. Our approach extends to 
medium-size nonlinear problems, due to automatic detection of the correlation sparsity pat¬ 
tern of input variables and roundoff error variables. Experimental results indicate that the 
optimization algorithm implemented in our Real2Float software package can often pro¬ 
duce bounds of quality similar to the ones provided by the competitive solvers Rosa and 
FPTaylor, while saving a significant amount of CPU time. In addition, Real2Float pro¬ 
duces sums of squares certificates which guarantee the correctness of these upper bounds 
and can can be efficiently verified inside the COQ proof assistant. 

This work yields several directions for further research investigation. First, we intend to 
increase the size of graspable instances by exploiting the SDP relaxations spe cifically tuned 
to th e case when a program implements the sum of many rational functions Bugarin et al. 
Symmetry patterns of certai n program sub-classes could be tackled with the SUP 

We could also provide roundoff error bounds for more 


2015 
hierarchies from 


Riener et al. 2013 

general programs, involving either finite or infinite conditional lo ops and additional compar- 

. A preliminary mandatory 
relaxations. Another inter- 
to derive sequences 


isons with the results described in Darulova and Kuncak 2016 


step is to be able to generate inductive invariants with SDP 


esting direction would be to apply the method used in Lasserre 2011 


of lower roundoff error bounds together with SDP-based certificates. On the formal proof 
side, we could benefit from floating-point/interval arithmetic libraries available inside COQ, 
first to improve the efficiency of the formal polynomial checker, currently relying on exact 
arithmetic, then to extend the formal verification to non-polynomial programs. The method 
implemented in FPTaylor happens to be efficient and precise to analyze various programs 
and it would be interesting to design a procedure combining FPTaylor with our tool on 
specific subsets of input constraints. Long-term research perspectives include theoretical 
study of why/when SOS performs better as well as satisfactory complexity results about 
SOS certificates (w.r.t. size of polynomials). Finally, we plan to combine this optimization 


framework with the procedure in 


to improve the automatic 


Gao and Constantinides 2015 
reordering of arithmetic expressions, allowing more efficient optimization of FPGA imple¬ 
mentations. 


NONLINEAR PROGRAM BENCHMARKS 


let box_rigidbody 1 x\ X2 x^, — [( — 15 , 15 ); ( — 15 , 15 ); (— 15 , 15 )]; ; 

let obj_rigidbodyl x\ X2 x^ — [{ — x\ X2 — 2 * X2 ^ x:^ — x\ — X'^^ Q)\\\ 

let box_rigidbody 2 x-i X2 tea — [( — 15 , 15 ); ( — 15 , 15 ); (— 15 , 15 )]; ; 

let ob j _r igidbody 2 X\X2X^ — 

[(2 * xi * 212 * + 3 * tea * tC 3 — 1 C 2 * * ^3 + 3 * tea * tcs — a: 2 , 0 )]; ; 

let box_keplerO xi X2 x^ X4 xq — [( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 )]; ; 
let obj_keplerO 2:1212^3 ^4 2152^6 = 

[(212 * X5 4 - X3 * xq — X2 * X3 — X5 * xq + Xi * {—Xi + 212 + 213 — 214 + 215 + 210), 0 )]; ; 

let box_keplerl 211 2:2 213 2:4 — [( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 )]; ; 

let obj_keplerl xi X2 x^ X4 — [{xi X4 ^ {—xi X2 x^ — X4) 

+212 * (211 — 212 + 213 + 214) + 213 * {xi + 212 — 213 + 214) 

— X2 * 213 * X4 — 211 * 213 — 211 * 212 — 214, 0)];; 

let box_kepler 2 211 2:2 213 2:4 215 xg = [( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 ); ( 4 , 6 . 36 )]; ; 

let obj_kepler2 xi X2 x^ X4 x^ xq — [{x\ * X4 ^ { — x\ X2 x:^ 

—X4 + 215 + 210) + 212 * 215 * {xi — 212 + 213 + 214 — X5 + 210) 
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+a:3 * rc6 * {^i X 2 — X 4 + — xq) — X 2 x^ * X 4 

— Xi * X^ X^ — Xi X 2 * Xq — X 4 * Xq * Xq, 0)]; ; 

let box_sineTaylor x — [(—1.57079632679, 1.57079632679)]; ; 
let obj_sineTaylor x — [(re — (re * re * rc)/6.0 
+ (rc*rc*rc*rc*rc)/120.0 

— (re *rc*rc*rc*rc*rc* a:)/5040.0, 0)]; ; 

let box_sine0rder3 2 =[( — 2,2)];; 

let obj _sine0rder3 2 = [(0.954929658551372 * 2 — 0.12900613773279798 * (2 * 2 * 2 ), 0)]; ; 
let box_sqroot y — [(0, 1)]; ; 

let obj_sqroot y = [(1.0 + 0.5*y — 0.125*y*y + 0.0625 ^ y ^ y y — 0.0390625 *y*-y*y*y, 0)];; 

let box_himmilbeau rci 3:2 = [( —5, 5); ( —5, 5)]; ; 
let obj _himmilbeau rci rc 2 = [( 

(rci * rci + rc2 — 11) * (rci * rci + rc2 — 11) + (rci + rc2 * rc2 — 7) * (rci + rc2 * rc2 — 7), 0)]; ; 
let box_dopplerl u v T = [(—100, 100); (20, 20e3); (—30, 50)]; ; 

let obj_dopplerl u v T = [(let ti=331.4 + 0.6*T in —ti*v/((ii+ii)*(ii+ii)),0)];; 

let box_doppler2 u v T = [(—125, 125); (15, 25e3); (—40, 60)]; ; 

let obj_doppler2 u v T = [(let ti=331.4 + 0.6*T in —ti*i;/((ii+ii)*(ii+ii)),0)];; 

let box_doppler3 u v T = [( — 30, 120); (320, 20300); ( — 50, 30)]; ; 

let obj_doppler3 u v T — [(let ti=331.4 + 0.6*T in —ti*i;/((ii+ii)*(ii+ii)),0)];; 

let box_verhulst rc = [(0.1,0.3)]; ; 

let obj_verhulst rc = [(4* x/(l + (rc/1.11)), 0)]; ; 

let box_carbonGas n = [(0.1, 0.5)]; ; 

let ob j _carbonGas v — [(let p — 3 . 5 e 7 in let a — 0.401 in 
let b — 42 . 7 e -6 in let t — 300 in let n = 1000 in 
(p -\- a * (n/v) * * 2 ) * ("U — n * 6) — 1 . 3806503 e -23 * n * i, 0 )]; ; 

let box_predPrey rc = [(0.1, 0.3)]; ; 

let obj_predPreyrc = [(4*rc*rc/(l + (rc/l.ll)**2), 0)];; 
let box_turbinel v w r = [(—4.5,—0.3); (0.4, 0.9); (3.8, 7.8)]; ; 

let obj_turbinel v w r=[(3 + 2/(r*r) — 0.125 *(3 — 2*n)*(tt;*ii;*r*r)/(l — t;) — 4.5,0)];; 

let box_turbine2 v w r = [(—4.5,—0.3); (0.4, 0.9); (3.8, 7.8)]; ; 

let obj_turbine2 v w r=[(6*'n — 0.5*1?* (w *it;*r*r)/(l — v) —2.5,0)];; 

let box_turbine3 v w r = [(—4.5,—0.3); (0.4, 0.9); (3.8, 7.8)]; ; 

let obj_turbine3 v w r=[(3—2/(r*r)—0.125*(l + 2*n)*(it;*ii;*r*r)/(l — i;) — 0.5,0)];; 

let box_jet rcirc 2 = [(—5, 5); ( —20, 5)]; ; 
let obj_jet rci rc 2 = [( 

rci + ((2*rci * ((3*rci *rci + 2*rc2 —3:i)/(rci *rci + 1)) 

*((3 * rci * rci + 2 * rc2 — rci)/(rci * rci + 1) — 3) 

+rci*rci*(4*((3*rci*rci+2*rc2—rci)/(rci*rci + l)) — 6)) 

*(rci * rci + 1) + 3 * rci * rci * ((3 * rci * rci + 2 * rc 2 — rci)/(rci * rci + 1)) 

+rci *rci *rci + rci + 3* ((3*rci *rci + 2*rc2 —rci)/(rci *rci + 1))),0)];; 

let box_f loudas2_6 rci rc 2 rcs rc4 3:5 rcg rcy rcs rcg rcio = 

[( 0 . 1 ); ( 0 , 1 ); ( 0 . 1 ); ( 0 . 1 ); ( 0 , 1 ); ( 0 . 1 ); ( 0 , 1 ); ( 0 . 1 ); ( 0 . 1 ); ( 0 , 1 )];; 

let cstr_f loudas2_6 xi X 2 xq X 4 xq xq xj xq xg xiq — [ 

—4 + 2*rci+6*3:2 + l*rc3 + 0*3:4 + 3*3?5+3*a?6 + 2*3:7 + 6*3?8 + 2*3?9 + 2* 3?io; 

22 — (6*3?!—5*3:2+8*3?3—3*3:4+0*3?5 + l*rc6+3*3:7 + 8*rc8+9*a;9 — 3* a;io); 

—6 — (5*a:i+6*3:2+5*3?3+3*rc4+8*a?5—8*a;6+9*3:7 + 2*rc8+0*3:9 — 9* a;io); 

— 23 — (9*3:i+5*3:2+0*a?3—9*rc4 + l*a? 5 —8*336+3*a?7 — 9*rc8~9* x9 — 3 * rcio)? 

— 12 — ( — 8 *3:i+7*3?2— 4 * 3:3 — 5 * 3 : 4 —9*a?5 + l*3:6 — 7*a?7 — 1*338+3*3:9 — 2* 2 : 10 )]; ; 
let obj_floudas2_6 xi X 2 X 3 X 4 xq xq xj xg xq xio — [( 

48 * rci + 42 * X 2 + 48 * 3:3 + 45 * 3:4 + 44 * 3:5 + 41 * 336 + 47 * xj 
+42 * 338 + 45 * 3:9 + 46 * rcio 

— 50 * (rci * 3:1 + 3:2 * 3:2 + 3:3 * 3:3 + 3:4 * 3:4 + 3:5 * 3:5 
+336 * Xq 4- xj * xj 4- xg * Xq + x9 * x9 4- xiQ * rcio), 0)];; 

let box_floudas3_3 3:1 3:2 3:3 3:4 3:5 srg = [(0,6);(0,6);(1,5);(0,6);(1,5);(0,10)];; 

— ( 3:2 — 2) * *2 — ( 3:3 — 1) * *2 — ( 3:4 — 4) * *2 — ( 3:5 — 1) * *2 — ( 3:3 — 4) * *2, 0)]; ; 


ACM Transactions on Mathematical Software, Vol. V, No. N, Article A, Publication date: January YYYY. 



A:28 


let cstr_f loudas3_3 xi X 2 X 4 ^ xq — 

[(x 3 — 3) * *2 + 3:4 — 4; { 2:5 — 3) * *2 + ^6 — 4; 

2 — £Ci + 3 * X 2 \ 2 + x\ — X 2 \ Q — x\ — X 2 \ x\ -\- X 2 — 2]; ; 

let obj_floudas3_3 x^ X 2 x:^ x^x^ xq — [{ — 2^ {xi — 2) ^ ^2 

— {x 2 — 2) * *2 — ( 3:3 — 1) * *2 — ( 3 J 4 — 4) * *2 

— ( 3:5 — 1) ^ ^2 — (3;e — 4) * *2, 0)];; 

let box_floudas3_4 xi3:2 3;3 — [(0,2);(0,2);(0,3)];; 
let cstr_f loudas3_4 3:1 3:;2 2:^3 = [ 

4 — 3^1 — 3^2 — 2 : 3 ; 6 — 3 * 3^2 — 3 : 3 ; 

— 0.75 + 2 * 3^1 — 2 * 3:3 + 4 * 3^1 * 3:1 — 4 * 3^1 * 3:2 

+4 * 3;i * 3^3 + 2 * 3^2 * 2:2 — 2 * 3:2 * 2:3 + 2 * 3^3 * 3 : 3 ]; ; 
let obj_floudas3_4 3 ?i 3 : 2 a :3 = [( — 2 * 3 :i+ 3:2 — 2 : 3 , 0 )];; 

let box_floudas4_6 3:1 3:2 — [(0,3); (0,4)]; ; 
let cstr_f loudas4_6 3:1 3:2 = [ 

2 * 3:1 * *4 — 8 * 3:1 * *3 + 8 * 3:1 * 3:1 — 3 : 2 ; 

4 * 3:1 * *4 — 32 * 3:1 * *3 + 88 * xi * 3:1 — 96 * 3:1 + 36 — 3 : 2 ];; 
let obj_floudas4_6 xx X 2 — [(— 2:1 — 2 : 2 , 0)]; ; 

let box_floudas 4 _ 73 :i 3:2 — [(0,2);(0, 3)];; 

let cstr_floudas4_7 2:1 3:2 = [— 2 * 3:1 **4 + 2 — 3 : 2 ];; 

let obj_floudas4_73:i3:2 = [(—12*3:i — 7* 3:2 + 2:2 * 3 : 2 , 0)];; 

let box_cavl0 3 : — [(0, 10)]; ; 

let obj_cavl0 3: = [( If {x x — x ^ 0) then 3:*0.1 else x*3: + 2, 0)];; 

let box_perin 3:y=[(l,7);( —2,7)];; 

let cstr_perin 3:2/=[3: — l;y + 2;3: — y;5—y — 2 :];; 

let obj_perin xy — [{ if (3:*x + y*?/<4) then y * x else 0,0)];; 

let box_logexp 3:=[( — 8 , 8 )];; 

let obj_logexp x — [(log(l + exp(3:)), 0)]; ; 

let box_sphere xryz = [(-10, 10); (0, 10); (-1.570796, 1.570796); (-3.14159265,3.14159265)]; ; 
let obj_sphere 3:ry2 = [(3: + r* sin(y) * cos( 2 ), 0)]; ; 

let box_hartman3 3:1 3:2 2:3 = [(0, 1); (0, 1); (0, 1)]; ; 
let obj_hartman3 2 : 13 : 22:3 = [( 

let el — 3.0 * ( 3:1 — 0.3689) * *2 + 10.0 * ( 2:2 — 0.117) * *2 + 30.0 * ( 3:3 — 0.2673) * *2 in 

let e2 — 0.1 * ( 3:1 — 0.4699) * *2 + 10.0 * ( 2:2 — 0.4387) * *2 + 35.0 * ( 2:3 — 0.747) * *2 in 

let e3 — 3.0 * ( 2:1 — 0.1091) * *2 + 10.0 * ( 3:2 — 0.8732) * *2 + 30.0 * ( 2:3 — 0.5547) * *2 in 

let e4 — 0.1 * ( 3:1 — 0.03815) * *2 + 10.0 * ( 2:2 — 0.5743) * *2 + 35.0 * ( 2:3 — 0.8828) * *2 in 

— (1.0 * exp(—el) + 1.2 * exp( —e2) + 3.0 * exp( —e3) + 3.2 * exp( —e4)), 0)]; ; 

let box_hartman 6 3:1 2:2 2:3 3:4 3:5 xq — [( 0 , 1 ); ( 0 , 1 ); ( 0 , 1 ); ( 0 , 1 ); ( 0 , 1 ); ( 0 , 1 )]; ; 
let obj_hartman 6 X 1 X 2 X 3 X 4 X 5 XQ — [( 

let el — 10.0 * ( 3:1 — 0.1312) * *2 + 3.0 * ( 2:2 — 0.1696) * *2 + 17. * ( 2:3 — 0.5569) * *2 + 3.5 * ( 3:4 — 0.0124) * *2 
+ 1.7 * ( 3:5 — 0.8283) * *2 + 8.0 * {xq — 0.5886) * *2 in 

let e2 — 0.05 * ( 3:1 — 0.2329) * *2 + 10 * ( 3:2 — 0.4135) * *2 + 17.0 * ( 2:3 — 0.8307) * *2 + 0.1 * ( 2:4 — 0.3736) * *2 

+8.0 * ( 3:5 — 0.1004) * *2 + 14.0 * ( 2:3 — 0.9991) * *2 in 

let e3 — 3.0 * ( 3:1 — 0.2348) * *2 + 3.5 * ( 3:2 — 0.1451) * *2 + 1.7 * ( 2:3 — 0.3522) * *2 + 10.0 * ( 3:4 — 0.2883) * *2 

+ 17.0 * ( 3:5 — 0.3047) * *2 + 8.0 * ( 2:3 — 0.665) * *2 in 

let e4 — 17.0 * ( 3:1 — 0.4047) * *2 + 8 * ( 2:2 — 0.8828) * *2 + 0.05 * ( 3:3 — 0.8732) * *2 + 10.0 * ( 2:4 — 0.5743) * *2 

+0.1 * ( 3:5 — 0.1091) * *2 + 14.0 * ( 2:3 — 0.0381) * *2 in 

— (1.0 * exp(—el) + 1.2 * exp( —e2) + 3.0 * exp( —e3) + 3.2 * exp( —e4)), 0)]; ; 
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